{
"cells": [
{
"cell_type": "markdown",
"id": "4f2ecd13",
"metadata": {},
"source": [
"# Populationsdynamik og Bæredygtighed\n",
"\n",
"Dette materiale er en tilpasset udgave af [Intermat 2.0 modul 3](https://intermat20.compute.dtu.dk/ct/Modul3_elever.html)."
]
},
{
"cell_type": "markdown",
"id": "3064cdce",
"metadata": {},
"source": [
"Det er helt centralt i forbindelse med politiske beslutninger for samfundet at kunne forudsige, hvordan en given population udvikler sig fremover. Det er selvsagt vigtigt at se på et lands befolkning og kunne planlægge dagsinstitutioner, skoler, ældrepleje etc., og vurdere den økonomiske bæredygtighed for den demografiske profil. Præcis samme overvejelser er meningsfulde for dyrearter. Det er for eksempel nødvendigt at kunne beskrive dynamikken for en fiskepopulation før man går i gang med at uddele fiskekvoter. Det er samme problemstillinger blot i meget forskellig kontekst.\n",
"\n",
"I denne notebook vil vi se på modellering af populationer ved brug af matricer og linear algebra. De grundlæggende ideer går tilbage til Patrick Leslie i 1945, som forelog at studere populationer ved at inddele i aldersklasser og trække på viden om klassernes individuelle fertilitets- og overlevelsesrater. Som eksempel så han på den brune rotte. Samlet giver det et diskret, dynamisk system, som kan analyseres med egenværdier og -vektorer.\n",
"\n",
"Nedenfor vil starte med små og simple modeller med udgangspunkt i kaninpopulationer og Fibonacci-tallene. Modellerne udvides, fertilitetsrater og overlevelsesrater introduceres. Endeligt når vi til en model for den danske befolknings udvikling med betingelser og rater baseret på data fra Danmarks statistik.\n",
"\n",
"Undervejs vil vi fremhæve væsentlige emner fra lineær algebra, matematisk modellering med Lesliematricer, egenværdier og -vektorer, iteration og konvergens, og vi vil tage nogle små sidespring for at belyse beregningselementer, herunder se på kontrolstrukturer og kompleksitet."
]
},
{
"cell_type": "markdown",
"id": "c7e1fbe3",
"metadata": {},
"source": [
"## Program\n",
"1. Opvarmning: Fibonacci-tallene genbesøgt\n",
"2. Getting real - fødsel og overlevelse\n",
"3. Den danske befolkning"
]
},
{
"cell_type": "markdown",
"id": "1eaf610c",
"metadata": {},
"source": [
"## Opsætning"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "25cab16c",
"metadata": {},
"outputs": [],
"source": [
"%matplotlib widget\n",
"from gym_cas import *\n",
"from spb import graphics, arrow_2d, geometry\n",
"from sympy import Circle, Point, re\n",
"from ipywidgets import FloatSlider, Layout\n",
"from numpy import array\n",
"from numpy.linalg import eig"
]
},
{
"cell_type": "markdown",
"id": "35b91f6a",
"metadata": {},
"source": [
"# 1. Opvarmning: Fibonacci-tallene genbesøgt\n",
"\n",
"Fibonacci var en italiensk matematiker i det 13. århundrede; han er også kendt som Leonardo Figlio di Bonacci. Gennem rejser til nord-Afrika stiftede han bekendtskab med det hindu-arabiske talsystem (tallene 0-9, decimaler), som han i sin bog Liber Abaci fra 1203 bragte til Europa. De nye tal afløste romertallene og bragte en revolution for handel og bankvæsenet; også dengang var matematik en kilde til magt og rigdom.\n",
"\n",
"Fibonacci havde et væsentligt fokus på algoritmer, ligninger og praktisk problem-løsning, og det er helt rimeligt at se ham som en pioner inden for computational thinking. Meget længe inden computerens fødsel. \n",
"\n",
"I Liber Abaci gennemfører Fibonacci et tankeeksperiment om en kaninpopulations udvikling i isolerede omgivelser. Udgangspunktet er et kaninpar og den matematiske model, at ethvert voksent par giver afkom i form af et nyt, ungt par hver måned. Unger bliver voksne på en måned, og der er ingen dødelighed. \n",
"\n",
"Vi vælger at tælle hunkaniner fremfor par. Vi starter med én voksen hunkanin, og i en kontekst af en populationsdynamik ser vi på en aldersopdelt populationsmodel, hvor der skelnes mellem to aldersklasser:\n",
"\n",
"- $\\boldsymbol{x}_k[0]$: antal **unge** hunkaniner efter $k$ måneder,\n",
"- $\\boldsymbol{x}_k[1]$: antal **voksne** hunkaniner efter $k$ måneder.\n",
"\n",
"Til $k=0$ har vi $\\boldsymbol{x}_0[0] = 0$ og $\\boldsymbol{x}_0[1] = 1.$ Dynamikken beskrives ved\n",
"\n",
"\\begin{align*}\n",
" \\boldsymbol{x}_{k+1}[0] &= \\boldsymbol{x}_{k}[1] &\\qquad \\text{(hver voksen kanin giver anledning til en unge i næste måned)},\\\\\n",
" \\boldsymbol{x}_{k+1}[1] &= \\boldsymbol{x}_{k}[0] + \\boldsymbol{x}_{k}[1] & \\qquad \\text{(hver unge giver anledning til en voksen i næste måned, og de nuværende voksne overlever)}.\n",
"\\end{align*}\n",
"\n",
"\n",
"Hvis vi introducerer matricen\n",
"\n",
"$$\n",
"A = \\begin{bmatrix}0 & 1\\\\[4pt]1 & 1\\end{bmatrix}\n",
"$$\n",
"så er populationen beskrevet ved det diskrete dynamiske system\n",
"\n",
"\\begin{align*}\n",
"\\boldsymbol{x}_{k+1} = A \\boldsymbol{x}_k, \\quad k = 0,1,2,\\ldots\\;\\;, \\quad \\boldsymbol{x}_0 = \\begin{bmatrix} 0 \\\\ 1 \\end{bmatrix}\n",
"\\end{align*}\n",
"\n",
"Dermed kan vi følge, hvordan antallet af unge og voksne udvikler sig måned for måned.\n",
"\n",
"Vi kan beskrive udviklingen ved at iterere matricen $A$ så\n",
"\n",
"$$\n",
"\\boldsymbol{x}_1 = A \\boldsymbol{x}_0, \\quad \n",
"\\boldsymbol{x}_2 = A \\boldsymbol{x}_1 = A^2 \\boldsymbol{x}_0, \\quad \\dots, \\quad\n",
"\\boldsymbol{x}_{k+1} = A \\boldsymbol{x}_k = A^{k+1} \\boldsymbol{x}_0.\n",
"$$\n",
"\n",
"Ser vi kun på det antal af voksne hunkaniner $F_k = \\boldsymbol{x}_k[1]$, genfinder vi rekursionen\n",
"\n",
"$$\n",
"F_{k+1} = F_k + F_{k-1}, \\quad k = 1,2,\\ldots\\;\\;,\n",
"$$\n",
"\n",
"som med begyndelsen $F_0 = F_1 = 1$ lige præcis giver **Fibonacci-tallene** $1,1,2,3,5,8,\\ldots$\n",
"\n",
"Matricen $A$ kaldes derfor ofte *Fibonacci-matricen*, da hver iteration svarer til én måned i modellen, og den voksen-bestanddel følger Fibonacci-følgen."
]
},
{
"cell_type": "markdown",
"id": "b78c0619",
"metadata": {},
"source": [
">Læg mærke til, hvordan første og anden række i $A$ beskriver henholdsvis næste generations unger og voksne. Fortolk ligeledes søjlerne i $A$ i en kontekst af kaninpopulationer."
]
},
{
"cell_type": "markdown",
"id": "2badd3bd",
"metadata": {},
"source": [
"## Iteration og for-løkker i Python\n",
"\n",
"Løkker (loops) er en såkaldt kontrolstruktur, som tillader os at løbe en gennem en sekvens af elementer og for hvert element udføre en given handling.\n",
"\n",
"Før vi går i gang med at implementere Fibonacci-matricen og iterere i Python, laver vi en simpel opvarmningsøvelse.\n",
"\n",
"\n",
"\n",
"Læs koden nedenfor, og svar på spørgsmålet **inden** du kører den i Python.\n",
"\n",
"> Hvad forventer du at der bliver printet? Kør koden og sammenlign med din forudsigelse."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "51a25ec9",
"metadata": {},
"outputs": [],
"source": [
"k = 1\n",
"for i in range(5):\n",
" k = k + i\n",
" print(\"i =\", i, \", k =\", k)"
]
},
{
"cell_type": "markdown",
"id": "c9ca43ec",
"metadata": {},
"source": [
"> Hvad sker der hvis vi ændrer `range(5)` til `range(1,6)`?\n",
"\n",
"> Færdiggør kodecellen nedenfor, så der benyttes et for-loop til at beregne og printe de første $10$ Fibonacci-tal. \\\n",
"> NB: I koden bruger vi variablerne `F0` og `F1` til løbende at beregne og overskrive de to seneste Fibonacci-tal. Det kan virke en smule forvirrende, fordi man kunne tro, at `F0` altid repræsenterer det første Fibonacci-tal og `F1` det næste. Her anvender vi dem dog kun som “pladsholdere” for de nyeste værdier i sekvensen - ikke som faste indeks."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1a3e1689",
"metadata": {},
"outputs": [],
"source": [
"F0 = \"INDSÆT KODE HER\"\n",
"F1 = \"INDSÆT KODE HER\"\n",
"\n",
"for i in range(\"INDSÆT KODE HER\"):\n",
" print(\"INDSÆT KODE HER\")\n",
" F0, F1 = \"INDSÆT KODE HER\" , \"INDSÆT KODE HER\""
]
},
{
"cell_type": "markdown",
"id": "1743c75f",
"metadata": {},
"source": [
"## Matrix-iteration"
]
},
{
"cell_type": "markdown",
"id": "fdee0c6f",
"metadata": {},
"source": [
"Vi har set, at kaninpopulationen kan beskrives med\n",
"\n",
"$$\n",
"\\boldsymbol{x}_{k+1} = A \\boldsymbol{x}_k = A^{k+1}\\boldsymbol{x}_0, \\quad \n",
"A = \\begin{bmatrix}0 & 1\\\\ 1 & 1\\end{bmatrix}.\n",
"$$\n",
"\n",
"hvor startpopulationen er\n",
"\n",
"$$\n",
"\\boldsymbol{x}_0 = \\begin{bmatrix}0\\\\ 1\\end{bmatrix} \\quad \\text{(0 unge, 1 voksen).}\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "86a97b8f",
"metadata": {},
"source": [
"> Færdiggør kodecellen nedenfor, så populationen for de første 10 måneder udregnes ved at gentage $\\boldsymbol{x}_{k+1} = A \\boldsymbol{x}_k$. Der printes antal unge, voksne og total for hvert måned. \\\n",
"> Du kan bruge `Matrix` til at definere $A$ og `vector` til at definere $x_0$."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f6c2936f",
"metadata": {},
"outputs": [],
"source": [
"A = \"INDSÆT KODE HER\"\n",
"\n",
"startpopulation = \"INDSÆT KODE HER\"\n",
"\n",
"antal_md = \"INDSÆT KODE HER\"\n",
"\n",
"print(\"Måned | Unge | Voksne | Total\")\n",
"population = startpopulation\n",
"for md in range(antal_md):\n",
" print(f\"{md:5d} | {int(population[0]):4d} | {int(population[1]):6d} | {int(sum(population)):3d}\")\n",
" # Opdater populationen\n",
" population = \"INDSÆT KODE HER\""
]
},
{
"cell_type": "markdown",
"id": "ea10808c",
"metadata": {},
"source": [
"I følgende kodecelle ønsker vi at vise udviklingen af voksne kaniner i et plot. Vi bruger listen `voksne_list` til løbende at gemme antallet i hver iteration. \n",
"> Færdiggør kodecellen nedenfor, så antallet af voksne kaniner afbildes for de første 10 måneder ved hjælp af matrix-iteration."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "091486e6",
"metadata": {},
"outputs": [],
"source": [
"A = \"INDSÆT KODE HER\"\n",
"\n",
"startpopulation = \"INDSÆT KODE HER\"\n",
"\n",
"antal_md = \"INDSÆT KODE HER\"\n",
"\n",
"# tome liste til antal voksne kaniner\n",
"voksne_list = []\n",
"\n",
"population = startpopulation.copy()\n",
"for md in range(antal_md):\n",
" # tilføj antal voksne kaniner til listen\n",
" voksne_list.append(population[1])\n",
"\n",
" # Opdater populationen\n",
" population = \"INDSÆT KODE HER\"\n",
"\n",
"\n",
"# Plot antal voksne kaniner\n",
"plot_points(\n",
" range(antal_md),\n",
" voksne_list,\n",
" xlabel=\"Måned\",\n",
" ylabel=\"Antal voksne kaniner\",\n",
" title=\"Udvikling af voksne kaniner (måned 0-9)\",\n",
")"
]
},
{
"cell_type": "markdown",
"id": "d90c5462",
"metadata": {},
"source": [
"> Udnyt at $\\boldsymbol{x}_{k} = A^{k}\\boldsymbol{x}_0$ til at finde populationen i måned 9 uden at iterere måned for måned. \\\n",
"> $A^{k}$ kan laves med `A**k`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "80fe8bb2",
"metadata": {},
"outputs": [],
"source": [
"# INDSÆT KODE HER"
]
},
{
"cell_type": "markdown",
"id": "d3f3d57a",
"metadata": {},
"source": [
"> Sammenlign dit resultat med opgaven, hvor du bestemte $\\boldsymbol{x}_{k+1}$ ved at lave gentagne $\\boldsymbol{x}_{k+1}=A\\boldsymbol{x}_k$."
]
},
{
"cell_type": "markdown",
"id": "659ce038",
"metadata": {},
"source": [
"## Egenværdier, singulærværdier og asymptotisk udvikling"
]
},
{
"cell_type": "markdown",
"id": "c2575a85",
"metadata": {},
"source": [
"En egenværdi $\\lambda \\in \\mathbb{R}$ og tilhørende egenvektor $\\boldsymbol{v} \\neq 0$ opfylder egenværdiproblemet\n",
"\n",
"$$\n",
"A \\boldsymbol{v} = \\lambda \\boldsymbol{v}.\n",
"$$\n",
"\n",
"Det betyder, at når vi anvender matricen $A$ på vektoren $\\boldsymbol{v}$, så bliver resultatet kun en skalering af $\\boldsymbol{v}$ med $\\lambda.$\n",
"\n",
"Geometrisk er $A v$ parallel med $\\boldsymbol{v}$, men længden ændres med faktor $\\lambda$.\n",
"\n",
"En egenværdi $\\lambda$ kan findes med determinanten og opfylder ligningen\n",
"\n",
"$$\n",
"\\det(A - \\lambda I) = 0.\n",
"$$\n",
"\n",
"For Fibonacci-matricen $A$ betyder det, at\n",
"\n",
"$$\n",
"\\det\\begin{pmatrix}-\\lambda & 1\\\\ 1 & 1-\\lambda\\end{pmatrix} = -\\lambda(1-\\lambda) - 1 = 0.\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "c8d2a0da",
"metadata": {},
"source": [
"> Løs andengradsligningen for $\\lambda$. \n",
"\n",
"> Indse, at der findes en positiv løsning, $\\lambda = \\lambda_1,$ og en negativ løsning $\\lambda = \\lambda_2,$ og at $\\lambda_1 > |\\lambda_2|.$ Indsæt værdierne i kodecellen nedenfor:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0b557019",
"metadata": {},
"outputs": [],
"source": [
"lambda1 = \"INDSÆT KODE HER\"\n",
"lambda2 = \"INDSÆT KODE HER\""
]
},
{
"cell_type": "markdown",
"id": "cd69927d",
"metadata": {},
"source": [
"Vi går nu til en geometrisk fortolkning.\n",
"\n",
"Nedenfor finder du den interaktive funktion `enhedsvektor_afbildning()`. Det ser omfattende ud, men du behøver ikke forstå hvert trin."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b56498df",
"metadata": {},
"outputs": [],
"source": [
"# Definer enhedsvektor som funktion af vinkel t\n",
"def xv(t):\n",
" return vector(cos(t), sin(t))\n",
"\n",
"\n",
"# Definer billedvektor y = A * x(t)\n",
"def yv(A, t):\n",
" return A * xv(t)\n",
"\n",
"\n",
"def enhedsvektor_afbildning(A):\n",
" \"\"\"\n",
" Interaktiv visualisering af enhedsvektorer x og deres billede y = Ax\n",
" for en 2x2-matrix A.\n",
" \"\"\"\n",
"\n",
" # Vinkler t fra 0 til 2pi\n",
" t_max = int(2 * 3.14 * 100)\n",
" ts = [t / t_max for t in range(t_max)]\n",
"\n",
" # Maksimal værdi for x/y for at skabe passende aksegrænser\n",
" max_val = max([max([*abs(xv(tt)), *abs(yv(A, tt))]) for tt in ts])\n",
"\n",
" params = {\n",
" t: FloatSlider(\n",
" value=0,\n",
" min=0,\n",
" max=2 * pi,\n",
" step=0.01,\n",
" description=\"t\",\n",
" layout=Layout(width=\"80%\"),\n",
" )\n",
" }\n",
" circle_plot = geometry(Circle(Point(0, 0), 1), fill=False, show_in_legend=False)\n",
" vx_plot = arrow_2d((0, 0), xv(t), \"\\\\vec{x}\", params=params)\n",
" vy_plot = arrow_2d((0, 0), yv(A, t), \"A\\\\vec{x}\", params=params)\n",
"\n",
" return graphics(\n",
" circle_plot,\n",
" vx_plot,\n",
" vy_plot,\n",
" aspect=\"equal\",\n",
" xlim=(-max_val, max_val),\n",
" ylim=(-max_val, max_val),\n",
" )"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "897b1d40",
"metadata": {},
"outputs": [],
"source": [
"# eksempel\n",
"A = Matrix([[0., 1],\n",
" [1, 1]])\n",
"enhedsvektor_afbildning(A)"
]
},
{
"cell_type": "markdown",
"id": "764d9522",
"metadata": {},
"source": [
"To vektorer $\\boldsymbol{v}_1$ og $\\boldsymbol{v}_2$ er lineært uafhængige, hvis\n",
"\n",
"$$\n",
"\\alpha \\boldsymbol{v}_1 + \\beta \\boldsymbol{v}_2 = \\boldsymbol{0}\n",
"\\ \\Rightarrow\\ \\alpha = 0 \\ \\text{og} \\ \\beta = 0.\n",
"$$\n",
"\n",
"> Aflæs ved brug af `enhedsvektor_afbildning()` to (lineært) uafhængige vektorer $\\boldsymbol{v}=\\boldsymbol{v}_1$ og $\\boldsymbol{v}=\\boldsymbol{v}_2$, for hvilke $A \\boldsymbol{v}$ er parallel med $\\boldsymbol{v}$. Brug evt. funktionen `xv(t)` efter du har fundet en passende værdi for $t$.\n",
"\n",
"Disse er egenvektorerne. Bemærk, at hvis $\\boldsymbol{v}$ er en egenvektor, så er $c\\boldsymbol{v}$ det også for ethvert tal $c\\in\\mathbb{R}$. Specielt er $-\\boldsymbol{v}$ også en egenvektor med samme egenværdi. \n",
"\n",
"> Aflæs tilsvarende de tilhørende egenværdier $\\lambda=\\lambda_1$ og $\\lambda=\\lambda_2.$\n",
"\n",
"> Definer dem i nedenstående kodecelle."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3afcd55a",
"metadata": {},
"outputs": [],
"source": [
"v1 = \"INDSÆT VÆRDIER HER\"\n",
"lamb1 = \"INDSÆT VÆRDIER HER\"\n",
"\n",
"v2 = \"INDSÆT VÆRDIER HER\"\n",
"lamb2 = \"INDSÆT VÆRDIER HER\""
]
},
{
"cell_type": "markdown",
"id": "69f8d248",
"metadata": {},
"source": [
">Sammenlign de aflæste egenværdier med de eksakte værdier $\\lambda_1,\\lambda_2$ fundet ovenfor."
]
},
{
"cell_type": "markdown",
"id": "5b8eb07a",
"metadata": {},
"source": [
"Egenvektorer hørende til en egenværdi $\\lambda$ kan findes ved at løse ligningen\n",
"\n",
"$$\n",
"(A-\\lambda I )\\boldsymbol{v} = 0.\n",
"$$\n",
"\n",
"I nedenstående kode bruger vi vores indsigt fra SVD til at finde en enhedsvektor $\\boldsymbol{w}_1,$ som (omtrent) er en egenvektor."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9a6fb6dd",
"metadata": {},
"outputs": [],
"source": [
"M = A - lambda1 * Matrix.eye(A.shape[0])\n",
"\n",
"# Beregn SVD af M\n",
"U, S, V = M.singular_value_decomposition()\n",
"\n",
"# Vi husker, at de singulære værdier i S er aftagende, så det sidste må afspejle nulrummet for M\n",
"w1 = V[:,1]"
]
},
{
"cell_type": "markdown",
"id": "c682b2ff",
"metadata": {},
"source": [
">Sammenlign den aflæste $\\boldsymbol{v}_1$ med fundne vektor $\\boldsymbol{w}_1$."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f8ce49d9",
"metadata": {},
"outputs": [],
"source": [
"display(w1,v1,S)"
]
},
{
"cell_type": "markdown",
"id": "b724202d",
"metadata": {},
"source": [
"I Sympy kan egenværdierne og deres tilhørende egenvektorer bestemmes med `A.eigenvects()` eller `eig(A)` fra Numpy."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5e22c16e",
"metadata": {},
"outputs": [],
"source": [
"(lamb2, _ , (v2,)), (lamb1, _ , (v1,)) = A.eigenvects()\n",
"display(lamb1, lamb2, v1, v2)"
]
},
{
"cell_type": "markdown",
"id": "def79adf",
"metadata": {},
"source": [
"> Sammenlign egenværdierne med de teoretiske værdier, som du tidligere har bestemt."
]
},
{
"cell_type": "markdown",
"id": "b2ce4324",
"metadata": {},
"source": [
"## Asymptotisk udvikling"
]
},
{
"cell_type": "markdown",
"id": "a11ce10b",
"metadata": {},
"source": [
"Vi kan skrive startvektoren som en linearkombination af egenvektorerne $\\boldsymbol{v_1}$ og $\\boldsymbol{v_2}$, så lad\n",
"\n",
"$$\n",
"\\boldsymbol{x}_0 = c_1 \\boldsymbol{v}_1 + c_2 \\boldsymbol{v}_2.\n",
"$$\n",
"\n",
"Når vi anvender matricen $A$ på $\\boldsymbol{x}_0$, får vi pga. linearitet, at\n",
"\n",
"$$\n",
"\\boldsymbol{x}_1 = A \\boldsymbol{x}_0 = A(c_1 \\boldsymbol{v}_1 + c_2 \\boldsymbol{v}_2) = c_1 A \\boldsymbol{v}_1 + c_2 A \\boldsymbol{v}_2 = c_1 \\lambda_1 \\boldsymbol{v}_1 + c_2 \\lambda_2 \\boldsymbol{v}_2.\n",
"$$\n",
"\n",
"Anvendes $A$ igen, fås\n",
"\n",
"$$\n",
"\\boldsymbol{x}_2 = A \\boldsymbol{x}_1 = A(c_1 \\lambda_1 \\boldsymbol{v}_1 + c_2 \\lambda_2 \\boldsymbol{v}_2)\n",
"= c_1 \\lambda_1 A \\boldsymbol{v}_1 + c_2 \\lambda_2 A \\boldsymbol{v}_2\n",
"= c_1 \\lambda_1^2 \\boldsymbol{v}_1 + c_2 \\lambda_2^2 \\boldsymbol{v}_2.\n",
"$$\n",
"\n",
"Gentager vi processen, får vi for den $k.$ iteration, at\n",
"\n",
"$$\n",
"\\boldsymbol{x}_k = c_1 \\lambda_1^k \\boldsymbol{v}_1 + c_2 \\lambda_2^k \\boldsymbol{v}_2.\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "26c302d1",
"metadata": {},
"source": [
"I det følgende ønsker vi at bruge egenvektor-dekompositionen til at beregne $\\boldsymbol{x}_k$ for $k=0,\\dots,K$ via\n",
"\n",
"$$\n",
"\\boldsymbol{x}_k = c_1 \\lambda_1^k \\boldsymbol{v}_1 + c_2 \\lambda_2^k \\boldsymbol{v}_2.\n",
"$$\n",
"\n",
"Først kræver det, at vi bestemmer konstanterne $c_1$ og $c_2$ for $\\boldsymbol{x}_0$. \n",
"\n",
"Brug, at vektoren \n",
"\n",
"$$\n",
"\\begin{bmatrix}c_1\\\\c_2\\end{bmatrix}\n",
"$$ \n",
"\n",
"løser ligningssystemet\n",
"\n",
"$$\n",
"\\begin{bmatrix}\\boldsymbol{v}_1&\\boldsymbol{v}_2\\end{bmatrix}\\begin{bmatrix}c_1\\\\c_2\\end{bmatrix}=\\boldsymbol{x}_0\n",
"$$\n",
"\n",
"Vink: Løsningen til et lineært ligningssystem $A\\boldsymbol{x}=\\boldsymbol{b}$ kan bestemmes ved `x=A.solve(b)`.
\n",
"Vink: Søljevektorer $\\boldsymbol{x}$ og $\\boldsymbol{y}$ kan samles til en matrix $M=\\begin{bmatrix}\\boldsymbol{x}&\\boldsymbol{y}\\end{bmatrix}$ ved `M=x.row_join(y)`.\n",
"\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "22c039e8",
"metadata": {},
"outputs": [],
"source": [
"startpopulation = vector(0., 1.)\n",
"# Saml systemmatrix\n",
"V = \"INDSÆT KODE HER\"\n",
"\n",
"# Løs ligningssystemet\n",
"c = \"INDSÆT KODE HER\"\n",
"\n",
"c1 = \"INDSÆT KODE HER\"\n",
"c2 = \"INDSÆT KODE HER\""
]
},
{
"cell_type": "markdown",
"id": "d2e20427",
"metadata": {},
"source": [
"Følgende celle genererer et plot for antal unge og voksne kaniner for $k=0,\\ldots,K$.\n",
"\n",
">Færdiggør cellen så $\\boldsymbol{x}_k$ beregnes via egenværdidekompositionen i hver iteration. \\\n",
">*Vink*: Man kan beregne $a^b$ i Python ved `a**b`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "13f57daf",
"metadata": {},
"outputs": [],
"source": [
"A = Matrix([[0.0, 1.0], [1.0, 1.0]])\n",
"\n",
"K = \"INDSÆT KODE HER\" # vælg fx 10 (beregn k=0..K)\n",
"\n",
"# Beregn x_k vha formel\n",
"unge = []\n",
"voksne = []\n",
"for k in range(K + 1):\n",
" xk = \"INDSÆT KODE HER\"\n",
" unge.append(xk[0])\n",
" voksne.append(xk[1])\n",
"\n",
"# Plot komponenterne som funktion af k\n",
"p1 = plot_points(\n",
" range(K + 1),\n",
" unge,\n",
" label=\"x[0] (unge)\",\n",
" show=False,\n",
" xlabel=\"Måned (k)\",\n",
" ylabel=\"Antal kaniner\",\n",
" title=\"Udvikling af hver aldersklasse over tid\",\n",
")\n",
"p2 = plot_points(range(K + 1), voksne, label=\"x[1] (voksne)\", show=False)\n",
"(p1 + p2).show()"
]
},
{
"cell_type": "markdown",
"id": "8aaa8a39",
"metadata": {},
"source": [
"> Tjek at de værdier du får for $\\boldsymbol{x}_k$ er de samme, som dem du tidligere har bestemt."
]
},
{
"cell_type": "markdown",
"id": "ccfcb2d8",
"metadata": {},
"source": [
"Hvis vi dividerer begge sider af egenvektor-dekompositionen med $\\lambda_1^k$, fås\n",
"\n",
"$$\n",
"\\frac{\\boldsymbol{x}_k}{\\lambda_1^k} = c_1 \\boldsymbol{v}_1 + c_2 \\left(\\frac{\\lambda_2}{\\lambda_1}\\right)^k \\boldsymbol{v}_2.\n",
"$$\n",
"\n",
"Da $|\\lambda_2 / \\lambda_1| < 1$, vil det andet led gå mod $0$ for $k\\to\\infty$, og dermed fås\n",
"\n",
"$$\n",
"\\frac{\\boldsymbol{x}_k}{\\lambda_1^k} \\to c_1 \\boldsymbol{v}_1 \\quad \\text{for} \\; k\\to \\infty.\n",
"$$\n",
"\n",
"Det betyder, at vektoren $\\boldsymbol{x}_k$ vokser med omtrent $\\lambda_1^k$, og retningen af $\\boldsymbol{x}_k$ nærmer sig $\\boldsymbol{v}_1$. \n",
"\n",
"\n",
"Hvis vi definerer den demografiske profil som\n",
"\n",
"$$\n",
"\\boldsymbol{d}_k = \\frac{\\boldsymbol{x}_k}{\\|\\boldsymbol{x}_k\\|},\n",
"$$\n",
"\n",
"så vis med papir og blyant, at når $\\boldsymbol{v}_1$ er en enhedsvektor, så er\n",
"\n",
"$$\n",
"\\boldsymbol{d}_k \\to \\boldsymbol{v}_1 \\quad \\text{for } k \\to \\infty.\n",
"$$\n",
"
\n",
"\n",
"\n",
"\n",
"💡 Hint (klik for at åbne)
\n",
"\n",
"Brug udtrykket \n",
"\n",
"$$\n",
"\\boldsymbol{x}_k = c_1 \\lambda_1^k \\boldsymbol{v}_1\n",
"+ c_2 \\lambda_2^k \\boldsymbol{v}_2\n",
"$$\n",
"\n",
"og faktoriser $\\lambda_1^k$ ud af både tæller og nævner i \n",
"\n",
"$$\n",
"\\boldsymbol{d}_k = \\frac{\\boldsymbol{x}_k}{\\|\\boldsymbol{x}_k\\|}\n",
"$$\n",
"\n",
"Udnyt til sidst, at\n",
"\n",
"$$\n",
"\\left(\\frac{\\lambda_2}{\\lambda_1}\\right)^k \\to 0 \\quad \\text{når } k\\to\\infty.\n",
"$$\n",
"\n",
" "
]
},
{
"cell_type": "markdown",
"id": "f191a6e3",
"metadata": {},
"source": [
"> Vis numerisk, at $\\boldsymbol{d}_k$ går mod $\\boldsymbol{v}_1$ for store $k$ ved at færdiggøre koden nedenfor. \\\n",
"> *Vink*: Hvis grafen ikke går mod 0, så prøv at sammenligne med $-\\boldsymbol{v}_1$."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a8588730",
"metadata": {},
"outputs": [],
"source": [
"K = 'INDSÆT KODE HER'\n",
"\n",
"# Beregn d_k = x_k / ||x_k|| for k = 0..K vha. egenvektor-formlen\n",
"dk_list = []\n",
"afstand = []\n",
"for k in range(K + 1):\n",
" xk = 'INDSÆT KODE HER'\n",
" dk = 'INDSÆT KODE HER'\n",
" dk_list.append(dk)\n",
" afstand.append((dk - abs(v1)).norm())\n",
"\n",
"plot_points(\n",
" range(K + 1),\n",
" afstand,\n",
" yscale=\"log\",\n",
" xlabel=\"k\",\n",
" ylabel=\"$\\\\|d_k - \\\\hat v_1\\\\|$\",\n",
" title=\"Afstand mellem $d_k$ og $\\\\hat v_1$\",\n",
")"
]
},
{
"cell_type": "markdown",
"id": "b4748fcc",
"metadata": {},
"source": [
"# 2. Getting real - fødsel og overlevelse \n",
"\n",
"Indtil nu har vi arbejdet med Fibonacci-matricen, hvor alle ikke-nul elementer var $1$ og populationen fulgte Fibonacci-følgen. Elementerne i matricen kan fortolkes som overlevelsesrater og fødselsrater. Vi generaliserer nu modellen ved at indføre.\n",
"\n",
"- $f_{2}$: fødselsrate - forventet antal unger per voksen per md. (sandsynlighed for afkom),\n",
"- $s_{1}$: overlevelsesrate fra ung til voksen (sandsynlighed for at en ung overlever til at blive voksen næste md), \n",
"- $s_{2}$: overlevelsesrate for voksne (sandsynlighed for at en voksen overlever til næste md).\n",
"\n",
"Den nye transitionsmatrix bliver:\n",
"\n",
"$$\n",
"A = \\begin{bmatrix}\n",
"0 & f_{2} \\\\\n",
"s_{1} & s_{2}\n",
"\\end{bmatrix}.\n",
"$$\n",
"\n",
"> Vis ved at se på den karakteristiske ligning $\\det(A-\\lambda I) = 0$ med papir og blyant, at hvis $f_2, s_1, s_2$ alle er positive, da har matricen to reelle egenværdier, $\\lambda_1>0$ og $\\lambda_2<0$, og at $\\lambda_1 > |\\lambda_2|.$"
]
},
{
"cell_type": "markdown",
"id": "f813d1fb",
"metadata": {},
"source": [
"\n",
"💡 Hint (klik for at åbne)
\n",
"\n",
"Den karakteristiske ligning kommer fra\n",
"\n",
"$$\n",
"\\det(A - \\lambda I) = 0.\n",
"$$\n",
"\n",
"For matricen\n",
"\n",
"$$\n",
"A = \\begin{bmatrix}\n",
"0 & f_2 \\\\\n",
"s_1 & s_2\n",
"\\end{bmatrix}\n",
"$$\n",
"\n",
"får du\n",
"\n",
"$$\n",
"\\det \\begin{bmatrix}\n",
"-\\lambda & f_2 \\\\\n",
"s_1 & s_2 - \\lambda\n",
"\\end{bmatrix}\n",
"= (-\\lambda)(s_2-\\lambda) - f_2 s_1 = 0.\n",
"$$\n",
"\n",
"Omskriv, så du får en andengradsligning i $\\lambda$ - og brug fortegnene på $f_2, s_1, s_2$ til at argumentere for fortegnene på rødderne samt hvilken der har størst absolutværdi.\n",
"\n",
" "
]
},
{
"cell_type": "markdown",
"id": "f988e221",
"metadata": {},
"source": [
"I følgende opgave arbejder vi ud fra samme fremgangsmåde som i første opgave, hvor vi blot nu undersøger en anden transitionsmatrix. Du kan genbruge størstedelen af koden, som vi brugte i første opgave.\n",
"\n",
"Opgaven går derfor ud på følgende:\n",
"\n",
"> 1. Vælg nogle værdier for $f_{2} \\in (0.2, 0,7), s_{1}, s_{2} \\in (0.7,1)$ og opbyg $A$. Overvej, hvad gode biologiske værdier kunne være (brug evt. en chatbot til at hjælpe dig). Husk at definere $A$ i en kodecelle.\n",
">\n",
"> 2. Vælg en startpopulation. Husk at definere den i en kodecelle.\n",
">\n",
"> 3. Brug `enhedsvektor_afbildning(A)`, `A.eigenvects()` eller `eig(A)` til at finde egenvektorerne $\\boldsymbol{v}_1$ og $\\boldsymbol{v}_2$ og egenværdier $\\lambda_1$ og $\\lambda_2$ som tidligere.\n",
">\n",
"> 4. Brug egenvektor-dekompositionen til at bestemme $x_k$ for $k=0,\\dots,K$ (brug formlen $\\boldsymbol{x}_k = c_1\\lambda_1^k \\boldsymbol{v}_1 + c_2\\lambda_2^k \\boldsymbol{v}_2$) og plot graferne for både $\\boldsymbol{x}_k[0]$ og $\\boldsymbol{x}_k[1]$.\n",
">\n",
"> 5. Vis numerisk, at $\\boldsymbol{d}_k$ stadig går mod $\\boldsymbol{v}_1$ for store $k$."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b6c35b85",
"metadata": {},
"outputs": [],
"source": [
"'INDSÆT KODE HER'"
]
},
{
"cell_type": "markdown",
"id": "77291523",
"metadata": {},
"source": [
"Når du har fået ovenstående til at fungere, kan du afprøve forskellige værdier - både for transitionsmatricen og for startpopulationen."
]
},
{
"cell_type": "markdown",
"id": "6397b4a0",
"metadata": {},
"source": [
"# 3. Den danske befolkning\n",
"\n",
"Vi udvider nu vores aldersstrukturerede modeller til en større population. Vi vælger at betragte hele den danske befolkning. \n",
"\n",
"- Populationen beskrives ved en vektor $\\boldsymbol{x} \\in \\mathbb{R}^{100}$, hvor hvert element $\\boldsymbol{x}[i]$ svarer til antallet af personer i alderen $i$. Vi betrager altså kun mennesker i alderen 0-99 år.\n",
"- Dynamikken styres af en **Leslie-matrix** $A \\in \\mathbb{R}^{100 \\times 100}$:\n",
"\n",
"$$\n",
"\\boldsymbol{x}_{k+1} = A \\boldsymbol{x}_k + \\boldsymbol{b}\n",
"$$\n",
"\n",
"Transitionsmatricen $A$ (der som nævnt er en Leslie-matrix) har følgende form:\n",
"\n",
"- Første række af $A$ indeholder **fødselsrater** $f_{15}, \\dots, f_{51} > 0$, altså aldersklasser hvor folk typisk får børn. \n",
"- Underdiagonal indeholder **overlevelsesrater** $s_i > 0$, som angiver sandsynligheden for at personer i aldersklasse $i$ overlever til aldersklasse $i+1$. \n",
"- Alle andre elementer i $A$ er nul. \n",
"\n",
"Dvs. $A$ har formen\n",
"\n",
"$$\n",
"A =\n",
"\\begin{bmatrix}\n",
"f_1 & f_2 & f_3 & \\cdots & f_n \\\\\n",
"s_1 & 0 & 0 & \\cdots & 0 \\\\\n",
"0 & s_2 & 0 & \\cdots & 0 \\\\\n",
"\\vdots & & \\ddots & \\ddots & \\vdots \\\\\n",
"0 & 0 & \\cdots & s_{n-1} & 0\n",
"\\end{bmatrix}.\n",
"$$\n",
"\n",
"$\\boldsymbol{b}$ kan bruges til at tilføje migration eller andre eksterne tilførsler."
]
},
{
"cell_type": "markdown",
"id": "454ffee3",
"metadata": {},
"source": [
"Vi ønsker at bruge den nuværende aldersfordeling af den danske befolkning som $\\boldsymbol{x}_0$.\n",
"\n",
"Brug filen `populationsdata.csv` eller hent nyeste aldersfordeling fra Danmarks Statistik ved at følge instruktionerne nedenfor:\n",
"\n",
"> 1. Åbn følgende link til [Danmarks Statistik](https://www.statistikbanken.dk/statbank5a/selectvarval/define.asp?PLanguage=0&subword=tabsel&MainTable=FOLK1AM&PXSId=239853&tablestyle=&ST=SD&buttons=0).\n",
"> 2. Vælg **Markér alle** under kategorien **ALDER**.\n",
"> 3. Klik på **VIS TABEL**.\n",
"> 4. Skift fra **Excel (\\*.xlsx)** til **Matrix (*.csv)** og klik på den blå pil lige ved siden af for at downloade CSV-filen.\n",
"> 5. Læg CSV-filen i samme directory som denne Jupyter-notebook.\n",
"> 6. Rediger `populationsdata.csv` til filnavnet på CSV-filen i nedenstående kodecelle."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ceeb184e",
"metadata": {},
"outputs": [],
"source": [
"# Indlæs CSV\n",
"with open(\"populationsdata.csv\", \"r\") as file:\n",
" file.readline()\n",
" aldre = []\n",
" startpopulation = []\n",
" for line in file:\n",
" alder, befolkning, dødstal, fertilitet, indvandring, udvandring = line.split(\n",
" \",\"\n",
" )\n",
" aldre.append(int(alder))\n",
" startpopulation.append(int(befolkning))"
]
},
{
"cell_type": "markdown",
"id": "1a53cabc",
"metadata": {},
"source": [
"Kør kodecellen nedenfor for at visualisere startpopulationen:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "82486842",
"metadata": {},
"outputs": [],
"source": [
"# Antal aldersklasser\n",
"k = len(startpopulation)\n",
"\n",
"# Plot startpopulation\n",
"plot_list(\n",
" aldre,\n",
" startpopulation,\n",
" xlabel=\"Alder\",\n",
" ylabel=\"Antal personer\",\n",
" title=\"Aldersfordeling i startpopulationen (0-99 år)\",\n",
")"
]
},
{
"cell_type": "markdown",
"id": "f61af573",
"metadata": {},
"source": [
"Vi ønsker at undersøge, hvordan befolkningen udvikler sig over tid ved hjælp af Leslie-matricen. Vi starter med en forenklet model, hvor fødselsrater og overlevelsesrater antages at være uniforme. \n",
"\n",
"> Definer Leslie-matricen. Vælg selv værdier for: \n",
"> - Uniform fødselsrate for aldersklasser 15-51 år.\n",
"> - Uniform overlevelsesrate for alle aldersklasser."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fedce3c7",
"metadata": {},
"outputs": [],
"source": [
"# Definér Leslie-matrix\n",
"k = 100\n",
"A = Matrix.zeros(k, k)\n",
"\n",
"# Fødselsrater\n",
"for i in range(15, 52):\n",
" A[0, i] = \"INDSÆT KODE HER\"\n",
"\n",
"# Overlevelsesrater\n",
"for i in range(k - 1):\n",
" A[i + 1, i] = \"INDSÆT KODE HER\""
]
},
{
"cell_type": "markdown",
"id": "26dd94d9",
"metadata": {},
"source": [
"I følgende kodecelle defineres en funktion `population_over_tid`, der simulerer en populations udvikling over tid givet en transitionmatrix, og en startpopulation.\n",
"\n",
"Funktionen generer også to figurer:\n",
"- total befolkning over tid.\n",
"- aldersfordelingen efter $K$ år.\n",
"\n",
">Kør cellen."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "efd7384f",
"metadata": {},
"outputs": [],
"source": [
"def population_over_tid(A, startpopulation, K, b = None, plot=True):\n",
" \"\"\"\n",
" Beregner population over K år givet transitionmatrix A og startpopulation.\n",
" Returnerer en matrix med population for hver aldersklasse over tid.\n",
" \"\"\"\n",
" k = A.shape[0]\n",
" X = Matrix.zeros(k,K+1)\n",
" X[:,0] = startpopulation\n",
"\n",
" # sætter b til nulvektor hvis ikke givet\n",
" if b is None:\n",
" b = Matrix.zeros(k,1)\n",
"\n",
" # iterationstep\n",
" for år in range(K):\n",
" X[:,år+1] = A * X[:,år] + b\n",
"\n",
" # generér plots hvis ønsket\n",
" if plot:\n",
" # plot total population over tid\n",
" total_population = Matrix.ones(1,k) * X\n",
" plot_points(range(K+1), total_population[:], xlabel='År', ylabel='Total befolkning', title='Udvikling af den danske befolkning')\n",
"\n",
" # plot aldersfordeling efter K år\n",
" plot_list(range(k), X[:,K][:], xlabel='Alder', ylabel='Antal personer', title=f'Aldersfordeling efter {K} år')\n",
"\n",
" return X"
]
},
{
"cell_type": "markdown",
"id": "189d722c",
"metadata": {},
"source": [
">Kør funktionen for din matrix $A$, og se plottene. Vælg selv antal år $K$."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "589f851e",
"metadata": {},
"outputs": [],
"source": [
"_ = population_over_tid(A, startpopulation, K=100)"
]
},
{
"cell_type": "markdown",
"id": "cc2fc109",
"metadata": {},
"source": [
"## Nettoreproduktionsrate\n",
"\n",
"**Nettoreproduktionsraten (NRR)** er det gennemsnitlige antal afkom, som en nyfødt individ forventes at producere i løbet af sin levetid. I praksis modelleres der som regel kun hunner, da de bestemmer den effektive reproduktion i populationen. Genkald formen af en Leslie-matrix\n",
"\n",
"$$\n",
"A =\n",
"\\begin{bmatrix}\n",
"f_1 & f_2 & f_3 & \\cdots & f_n \\\\\n",
"s_1 & 0 & 0 & \\cdots & 0 \\\\\n",
"0 & s_2 & 0 & \\cdots & 0 \\\\\n",
"\\vdots & & \\ddots & \\ddots & \\vdots \\\\\n",
"0 & 0 & \\cdots & s_{n-1} & 0\n",
"\\end{bmatrix}.\n",
"$$\n",
"\n",
"For en $n\\times n$ Leslie-matrix beregnes nettoreproduktionsraten $R_0$ som\n",
"\n",
"$$\n",
"R_0 = \\sum_{i=1}^{n} f_i \\, l_i ,\n",
"$$\n",
"\n",
"hvor $f_i$ er fertilitetsraten for aldersklasse $i$, og $l_i$ er sandsynligheden for at overleve til begyndelsen af aldersklasse $i$. Overlevelsessandsynligheden $l_i$ defineres rekursivt som\n",
"\n",
"$$\n",
"l_1 = 1, \\qquad\n",
"l_2 = s_1, \\qquad\n",
"l_3 = s_1 s_2, \\qquad \\dots, \\qquad\n",
"l_i = \\prod_{j=1}^{i-1} s_j ,\n",
"$$\n",
"\n",
"hvor $s_j$ er sandsynligheden for at overleve fra aldersklasse $j$ til $j+1$. Dermed kan $R_0$ også skrives eksplicit som\n",
"\n",
"$$\n",
"R_0 = f_1 + f_2 s_1 + f_3 s_1 s_2 + \\dots + f_n s_1 s_2 \\cdots s_{n-1} =\\sum_{i=1}^n f_i\\left( \\prod_{j=1}^{i-1}s_j \\right).\n",
"$$\n",
"\n",
"I følgende kodecelle definerer vi en funktion, der beregner nettoreproduktionsraten for en given Leslie-matrix.\n",
"\n",
">Færdiggør funktionen.\n",
"\n",
"\n",
"💡 Hint (klik for at åbne)
\n",
"\n",
"Tænk på, hvordan værdierne i en Leslie-matrix er organiseret:\n",
"\n",
"- Den første række indeholder **fertilitetsrater** for hver aldersklasse.\n",
"- De andre rækker indeholder **overlevelsesrater**, som viser sandsynligheden for at gå fra én aldersklasse til den næste.\n",
"\n",
"Overvej, hvordan du kan hente disse værdier i en Python-matrice med `A[i,j]` og bruge dem i løkkerne for `fi` og `sj`.\n",
"\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c22fe931",
"metadata": {},
"outputs": [],
"source": [
"def NRR(A):\n",
" \"\"\"\n",
" Beregner nettoreproduktionsraten R0 for Leslie-matrix A.\n",
" \"\"\"\n",
" n = A.shape[0]\n",
" R0 = 0.0\n",
"\n",
" # Beregn R0\n",
" for i in range(n):\n",
" # fødselsrate for aldersklasse i\n",
" fi = \"INDSÆT KODE HER\"\n",
"\n",
" # overlevelsesrate fra fødsel til aldersklasse i\n",
" li = 1.0\n",
" for j in range(i):\n",
" # overlevelsesrate for aldersklasse j til j+1\n",
" sj = \"INDSÆT KODE HER\"\n",
" li *= sj\n",
" R0 += fi * li\n",
"\n",
" return R0"
]
},
{
"cell_type": "markdown",
"id": "99a567e0",
"metadata": {},
"source": [
"Fortolkningen af $R_0$ er, at \n",
"- hvis $R_0>1$ forventes populationen af vokse.\n",
"- hvis $R_0=1$ er populationen stabil.\n",
"- hvis $R_0<1$ forventes populationen at falde.\n",
"\n",
"> Brug funktionen til at undersøge nettoreproduktionsraten. Svarer fortolkningen til hvad du kunne se af plottene?"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7e13e571",
"metadata": {},
"outputs": [],
"source": [
"\"INDSÆT KODE HER\""
]
},
{
"cell_type": "markdown",
"id": "09a3cd94",
"metadata": {},
"source": [
"Vi har de samme fortolkninger for den dominerende egenværdi, og om den er over/under 1."
]
},
{
"cell_type": "markdown",
"id": "8d0d6329",
"metadata": {},
"source": [
"I den følgende kodecelle definerer vi en funktion `dominating_eigenvalue`, der returnerer den dominerende egenværdi og den tilhørende egenvektor for en given matrix $A$."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2b6d3b35",
"metadata": {},
"outputs": [],
"source": [
"def dominating_eigenvalue(A):\n",
" \"\"\"\n",
" Beregner den dominerende egenværdi og tilhørende egenvektor for matrix A.\n",
" \"\"\"\n",
" eigenvalues, eigenvectors = eig(array(A).astype(\"float\"))\n",
"\n",
" # Finder indeks for den største egenværdi (absolut værdi)\n",
" idx = max(range(len(eigenvalues)), key=lambda i: eigenvalues[i])\n",
"\n",
" # største egenværdi-par\n",
" largest_eigenvalue = eigenvalues[idx]\n",
" largest_eigenvector = vector(*eigenvectors[:, idx])\n",
"\n",
" # Normaliserer til norm 1 og positiv sum\n",
" largest_eigenvector = largest_eigenvector / sum(largest_eigenvector)\n",
" largest_eigenvector = largest_eigenvector / largest_eigenvector.norm()\n",
"\n",
" # tager realdel\n",
" largest_eigenvalue = re(largest_eigenvalue)\n",
" largest_eigenvector = re(largest_eigenvector)\n",
"\n",
" return largest_eigenvalue, largest_eigenvector\n",
"\n",
"# Print results\n",
"largest_eigenvalue, largest_eigenvector = dominating_eigenvalue(A)\n",
"display(\"Largest eigenvalue:\", largest_eigenvalue)\n",
"display(\"Corresponding eigenvector:\", largest_eigenvector)"
]
},
{
"cell_type": "markdown",
"id": "858830d6",
"metadata": {},
"source": [
">Undersøg den dominerende egenværdi vha. `dominating_eigenvalue`.\n",
"\n",
">Svarer fortolkningerne til hinanden? Er $R_0=\\lambda_{max}$?"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3bd2f79e",
"metadata": {},
"outputs": [],
"source": [
"\"INDSÆT KODE HER\""
]
},
{
"cell_type": "markdown",
"id": "8f611758",
"metadata": {},
"source": [
">Prøv at gå tilbage til hvor Leslie-matricen blev defineret. Kan du finde værdier for fødsel- eller overlevelsesraterne, hvor nettoreproduktionsraten er tæt på 1? Vis plottene igen, og se den langvarige udvikling."
]
},
{
"cell_type": "markdown",
"id": "ba999242",
"metadata": {},
"source": [
"## Ny fødselsrate som parabel"
]
},
{
"cell_type": "markdown",
"id": "08377e06",
"metadata": {},
"source": [
"Vi vil nu eksperimentere med en ikke-uniform fødselsrate, som varierer med alderen. I stedet for at vælge en konstant rate for alle fertile aldersklasser, lader vi raten følge en parabel, der har top i alder 32 og nul i alderen 14 og 52.\n",
"\n",
"> Brug denne parabel til modellering af fødselsraten til at opdatere Leslie-matrixen i koden nedenfor, hvor du igen får\n",
"> - Et plot for total befolkning over tid.\n",
"> - Et plot, der viser, hvordan populationen fordeler sig på aldersklasser efter $K$ år.\n",
">\n",
"> Du kan selv eksperimentere med, hvad toppunktsværdien skal være.\n",
"\n",
"\n",
"💡 Hint (klik for at åbne)
\n",
"\n",
"Tænk på fødselsraten som en parabel, der er nul ved de yngste og ældste fertile aldersklasser og har et maksimum ved den alder, hvor fertiliteten er størst.\n",
"\n",
"$$\n",
"f(i) = f_{\\max} \\cdot \\frac{(i - i_\\text{min})(i_\\text{max} - i)}{(i_\\text{top} - i_\\text{min})(i_\\text{max} - i_\\text{top})}\n",
"$$\n",
"\n",
"hvor $i_\\text{min}$ og $i_\\text{max}$ er de yderste fertile aldre, $i_\\text{top}$ er toppen, og $f_{\\max}$ er den maksimale rate.\n",
"\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "963f7b3a",
"metadata": {},
"outputs": [],
"source": [
"# Antal aldersklasser\n",
"k = 100\n",
"\n",
"# Definér Leslie-matrix\n",
"A = Matrix.zeros(k, k)\n",
"\n",
"# Fødselsrater: udfyld parablen her\n",
"for i in range(15, 52):\n",
" # beregn A[0, i] som parabel med top i 32 og nul i 14 og 52\n",
" A[0, i] = \"INDSÆT KODE HER\"\n",
"\n",
"# Overlevelsesrater\n",
"for i in range(k - 1):\n",
" A[i + 1, i] = 0.95\n",
"\n",
"# Beregn population over K år\n",
"print(NRR(A))\n",
"X = population_over_tid(A, startpopulation, K=100)"
]
},
{
"cell_type": "markdown",
"id": "5f96dacf",
"metadata": {},
"source": [
"## Ny dødsrate som lineær funktion"
]
},
{
"cell_type": "markdown",
"id": "48bbc957",
"metadata": {},
"source": [
"Vi vil nu udvide vores model af den danske befolkning yderligere. \n",
"- Fødselsraten skal være en parabel, som vi lavede tidligere (top ved 32 år, nulpunkt ved 14 og 50 år). \n",
"- Nyt tillæg: Dødsraten skal være en lineær funktion med laveste dødsrate ved 0 år, stigende mod ældre aldersklasser."
]
},
{
"cell_type": "markdown",
"id": "edea3551",
"metadata": {},
"source": [
"> Indfør den nye dødsrate i koden nedenfor. Du vil igen få:\n",
"> - Et plot, der viser den samlede befolkning over tid.\n",
"> - Et plot, der viser aldersfordelingen af populationen efter $K$ år.\n",
">\n",
"> Du kan selv eksperimentere med dødsraten i år 0 og år 99.\n",
"\n",
"\n",
"\n",
"💡 Hint (klik for at åbne)
\n",
"\n",
"En lineær funktion kan skrives generelt som\n",
"\n",
"$$\n",
"y = y_\\text{max} - (y_\\text{max}-y_\\text{min}) \\frac{x}{x_\\text{max}}.\n",
"$$\n",
"\n",
"Formlen sikrer, at overlevelsesraten falder **lineært** fra høj værdi ved unge til lav værdi ved gamle. \n",
"Du kan bruge denne lineære formel inde i løkken for `A[i+1, i]`.\n",
"\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0d4ef8a8",
"metadata": {},
"outputs": [],
"source": [
"# Antal aldersklasser\n",
"k = 100\n",
"\n",
"# Definér Leslie-matrix\n",
"A = Matrix.zeros(k, k)\n",
"\n",
"# Fødselsrater\n",
"for i in range(15, 52):\n",
" A[0, i] = \"INDSÆT KODE HER\"\n",
"\n",
"# Lineær overlevelsesrate\n",
"for i in range(k - 1):\n",
" A[i + 1, i] = \"INDSÆT KODE HER\"\n",
"\n",
"# Beregn population over K år\n",
"print(NRR(A))\n",
"X = population_over_tid(A, startpopulation, K=100)"
]
},
{
"cell_type": "markdown",
"id": "15d6dfee",
"metadata": {},
"source": [
"## Beregningskompleksitet af matrix-iteration\n",
"\n",
"\n",
"Vi har tidligere i opgaven undersøgt populationen ved hjælp af Leslie-matricen og beregnet populationen over tid ved iterativt at anvende\n",
"\n",
"$$\n",
"\\boldsymbol{x}_{k+1} = A \\boldsymbol{x}_k = A (A (... (A \\boldsymbol{x}_0) ...)).\n",
"$$\n",
"\n",
"Vi har altså **ikke** beregnet $A^{k+1} \\boldsymbol{x}_0$ direkte, men gentaget matrix-vektor-multiplikationen $k+1$ gange.\n",
"\n",
"I denne opgave skal vi estimere, hvor mange grundlæggende beregninger (multiplikationer og additioner) der kræves for en $100 \\times 100$ Leslie-matrix.\n",
"\n",
"- Brug $k = 50$ \n",
"- Dimension: $A$ er $100 \\times 100$ \n",
"\n",
"Vi definerer først $k$ og dimensionen $n$."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f7fd4930",
"metadata": {},
"outputs": [],
"source": [
"k = 50\n",
"n = 100"
]
},
{
"cell_type": "markdown",
"id": "64caf7ff",
"metadata": {},
"source": [
">Overvej for matrix-vektor produktet:\n",
">- Hvor mange multiplikationer bruges?\n",
">- Hvor mange additioner bruges?\n",
">\n",
">*Vink*: For dot-produktet mellem to vektorer i $\\mathbb{R}^n$ bruges $n$ multiplikationer og $n-1$ additioner.\n",
"\n",
">Definer resultatet i følgende celle.\\\n",
">*Vink*: I Python skrives $x^a$ som `x**a`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1a313c74",
"metadata": {},
"outputs": [],
"source": [
"mv_mul = \"INDSÆT KODE HER\"\n",
"mv_add = \"INDSÆT KODE HER\"\n",
"mv_ops = mv_mul + mv_add"
]
},
{
"cell_type": "markdown",
"id": "77c6934b",
"metadata": {},
"source": [
">Overvej det samme for et matrix-matrix produkt og definer resultatet."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "75da657b",
"metadata": {},
"outputs": [],
"source": [
"mm_mul = \"INDSÆT KODE HER\"\n",
"mm_add = \"INDSÆT KODE HER\"\n",
"mm_ops = mm_mul + mm_add"
]
},
{
"cell_type": "markdown",
"id": "b4fddb89",
"metadata": {},
"source": [
">Hvor mange hhv. matrix-matrix og matrix-vektor produkter anvendes til at regne $(A^{k+1})\\boldsymbol{x}_0$ direkte? \n",
"\n",
">Brug dette til at finde det totale antal operationer i udregningen og definer resultatet."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "51d90814",
"metadata": {},
"outputs": [],
"source": [
"direkte_total = \"INDSÆT KODE HER\""
]
},
{
"cell_type": "markdown",
"id": "128df8d7",
"metadata": {},
"source": [
">Overvej på samme måde antallet af hhv. matrix-matrix og matrix-vektor produkter, der bruges til den iterative metode $\\boldsymbol{x}_{k+1}=A(A(...A(A\\boldsymbol{x}_0)))$.\n",
"\n",
">Definer det totale antal operationer."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a1e41416",
"metadata": {},
"outputs": [],
"source": [
"iterativ_total = \"INDSÆT KODE HER\""
]
},
{
"cell_type": "markdown",
"id": "ecc65bd9",
"metadata": {},
"source": [
"Følgende celle printer de to resultater og forholdet mellem dem.\n",
"\n",
">Kør cellen."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "42b43c0d",
"metadata": {},
"outputs": [],
"source": [
"print(\"Direkte metode:\\n\",direkte_total)\n",
"print(\"Iterativ metode:\\n\",iterativ_total)\n",
"print(\"Sammenligning:\\nDirekte/Iterativ operationer:\",direkte_total / iterativ_total)"
]
}
],
"metadata": {
"jupytext": {
"text_representation": {
"extension": ".md",
"format_name": "myst",
"format_version": 0.13,
"jupytext_version": "1.17.2"
}
},
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.0"
},
"source_map": [
12,
16,
20,
30,
37,
41,
59,
113,
117,
131,
136,
143,
150,
154,
169,
174,
187,
192,
219,
224,
226,
230,
234,
258,
264,
267,
273,
348,
353,
370,
376,
380,
390,
398,
402,
406,
410,
414,
418,
422,
450,
477,
487,
494,
518,
522,
575,
581,
603,
624,
656,
672,
674,
678,
712,
725,
733,
737,
752,
760,
772,
782,
821,
825,
827,
880,
902,
911,
913,
917,
921,
949,
955,
957,
961,
965,
986,
1004,
1008,
1014,
1035,
1052,
1072,
1075,
1086,
1090,
1094,
1098,
1104,
1106,
1112,
1114,
1120
]
},
"nbformat": 4,
"nbformat_minor": 5
}